new_results_analysis

library("reticulate")

## Warning: package 'reticulate' was built under R version 4.0.5

library("knitr")
library("Hmisc")

## Loading required package: lattice

## Loading required package: survival

## Loading required package: Formula

## Loading required package: ggplot2

## 
## Attaching package: 'Hmisc'

## The following objects are masked from 'package:base':
## 
##     format.pval, units

library("DescTools")

## 
## Attaching package: 'DescTools'

## The following objects are masked from 'package:Hmisc':
## 
##     %nin%, Label, Mean, Quantile

library("stringr")
library("egg")

## Loading required package: gridExtra

library("tidyverse")

## ── Attaching packages ─────────────────────────────────────── tidyverse 1.3.1 ──

## ✓ tibble  3.1.6     ✓ purrr   0.3.4
## ✓ tidyr   1.1.4     ✓ dplyr   1.0.7
## ✓ readr   2.1.1     ✓ forcats 0.5.1

## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## x dplyr::combine()   masks gridExtra::combine()
## x dplyr::filter()    masks stats::filter()
## x dplyr::lag()       masks stats::lag()
## x dplyr::src()       masks Hmisc::src()
## x dplyr::summarize() masks Hmisc::summarize()

# set plotting theme 
theme_set(theme_classic() + 
    theme(text = element_text(size = 24)))

# knitr chunk display options     
opts_chunk$set(comment = "",
               results = "hold",
               fig.show = "hold")

# suppress summarise() grouping warning 
options(dplyr.summarise.inform = F)

Load Data

use_condaenv("plinko")
pd = import("pandas")
df_data = pd$read_pickle("../../data/human_data/full_dataset_vision_corrected.xz")

Segment judgment, rt, and eye-data

df_data_judge = df_data %>% 
  select(participant, trial, response) %>% 
  unique()

df_data_rt = df_data %>% 
  group_by(participant, trial) %>% 
  summarise(rt = tail(t, n=1) - head(t, n=1)) %>% 
  mutate(log_rt = log(rt))

Compute Judgment Means

df_data_mean_judge_train = df_data_judge %>% 
  filter(participant %in% seq(1,15)) %>% 
  group_by(trial) %>% 
  summarise(hole1 = sum(response == 1)/n(),
            hole2 = sum(response == 2)/n(),
            hole3 = sum(response == 3)/n()) %>% 
  pivot_longer(c(hole1, hole2, hole3), 
               names_to = "hole",
               values_to = "human_mean")

Compute RT Means

df_data_mean_rt_train = df_data_rt %>% 
  filter(participant %in% seq(1,15)) %>% 
  mutate(log_rt = ifelse(rt != 0, log(rt), 0)) %>% 
  group_by(trial) %>% 
  summarise(mean_rt = mean(rt),
            mean_log_rt = mean(log_rt))

Parameter Search

compare_bandit_judgment = function(df_human_means, df_model_raw) {
  
  df_model_means = df_model_raw %>% 
    select(trial, run, judgment) %>% 
    group_by(trial) %>% 
    summarise(hole1 = sum(judgment == 0)/n(),
              hole2 = sum(judgment == 1)/n(),
              hole3 = sum(judgment == 2)/n()) %>% 
    pivot_longer(c(hole1, hole2, hole3),
                 names_to = "hole",
                 values_to = "model_mean")
  
  sq_err = left_join(df_human_means,
                     df_model_means,
                     by = c("trial", "hole")) %>% 
    summarise(sq_err = sum((human_mean - model_mean)^2)) %>% 
    pull(sq_err)
  
  return(sq_err)
  
}

compare_bandit_rt = function(df_human_means, df_model_raw) {
  df_model_means = df_model_raw %>% 
    select(trial, run, num_cols) %>% 
    mutate(time_measure = num_cols,
           log_time = ifelse(num_cols != 0, log(num_cols), num_cols)) %>% 
    group_by(trial) %>% 
    summarise(mean_time = mean(time_measure),
              mean_log_time = mean(log_time))
  
  df_compare = left_join(df_human_means,
                     df_model_means,
                     by = c("trial")) 
  
  
  model = lm(mean_log_rt ~ mean_log_time,
             data = df_compare)
  
  scaled_model_predictions = predict(model)
  
  sq_err = sum((df_compare$mean_log_rt - scaled_model_predictions)^2)
  
  return(sq_err)
}

compare_fixed_judgment = function(df_human_means, df_model_raw) {
  df_model = df_model_raw %>%
    select(trial, hole1, hole2, hole3) %>%
    pivot_longer(c(hole1, hole2, hole3),
                 names_to = "hole",
                 values_to = "model_score")
  
  
  sq_err = left_join(df_human_means,
                     df_model,
                     by = c("trial", "hole")) %>% 
    summarise(sq_err = sum((human_mean - model_score)^2)) %>% 
    pull(sq_err)
  
  return(sq_err)
}

compare_fixed_rt = function(df_human_means, df_model_raw) {
  df_model = df_model_raw %>% 
    select(trial, num_cols) %>% 
    mutate(log_time_measure = log(num_cols))
  
  df_compare = left_join(df_human_means,
                     df_model,
                     by = c("trial"))
  
  
  model = lm(mean_log_rt ~ log_time_measure,
             data = df_compare)
  
  scaled_model_predictions = predict(model)
  
  sq_err = sum((df_compare$mean_log_rt - scaled_model_predictions)^2)
  
  return(sq_err)
}

Load Bandit Heatmap Performance

path = "../python/model/model_performance/grid_regression_error/"
reg_files = list.files(path)

df_reg_perf = tibble()

for (file in reg_files) {
  if (grepl("bandit", file)) {
    temp = read.csv(paste(path, file, sep = ""))
    df_reg_perf = rbind(df_reg_perf, temp)
  }
}

df_reg_perf = df_reg_perf %>% 
  rename(thresholds = threshold,
         tradeoffs = tradeoff,
         sample_weights = sample_weight,
         look_dist = sq_err)

thresholds = c()
tradeoffs = c()
bandit_bws = c()
sample_weights = c()
bandit_judge_err = c()
bandit_rt_err = c()
bandit_look_dist = c()


path = "../python/model/model_performance/grid_judgment_rt/"
list_filenames = list.files(path)

for (file in list_filenames) {
  file_list = str_split(file, "_")[[1]]
  if ("bandit" %in% file_list) {
    
    threshold = as.numeric(file_list[match("threshold", file_list) + 1])
    tradeoff = as.numeric(file_list[match("tradeoff", file_list) + 1])
    bandit_bw = as.numeric(file_list[match("bw", file_list) + 1])
    sample_weight = as.numeric(file_list[match("weight", file_list) + 1])
    
    df_model_raw = read.csv(paste(path, file, sep = ""))
    
    judge_err = compare_bandit_judgment(df_data_mean_judge_train, df_model_raw)
    rt_err = compare_bandit_rt(df_data_mean_rt_train, df_model_raw)
    
    thresholds = append(thresholds, threshold)
    tradeoffs = append(tradeoffs, tradeoff)
    bandit_bws = append(bandit_bws, bandit_bw)
    sample_weights = append(sample_weights, sample_weight)
    bandit_judge_err = append(bandit_judge_err, judge_err)
    bandit_rt_err = append(bandit_rt_err, rt_err)
    
  }
  
  
}

df_bandit_performance = tibble(thresholds = thresholds,
                              tradeoffs = tradeoffs,
                              bws = bandit_bws,
                              sample_weights = sample_weights,
                              judge_err = bandit_judge_err,
                              rt_err = bandit_rt_err)

df_bandit_performance = df_bandit_performance %>% 
  left_join(df_reg_perf,
            by = c("thresholds",
                   "tradeoffs",
                   "bws",
                   "sample_weights"))

Bandit

df_bandit_performance = df_bandit_performance %>% 
  mutate(rank_judge = rank(judge_err),
         rank_rt = rank(rt_err),
         rank_look_dist = rank(look_dist),
         combined_ranks = (rank_judge + rank_rt + rank_look_dist)/3) %>% 
  arrange(combined_ranks)

df_to_show = df_bandit_performance %>%
  filter(judge_err < 38,
         rt_err < 10000)

g1 = ggplot(df_to_show, mapping = aes(judge_err, rt_err)) +
  geom_point() +
  geom_smooth(method = "lm") +
  xlab("Judgment Error") +
  ylab("RT Error")

g2 = ggplot(df_to_show, mapping = aes(judge_err, look_dist)) +
  geom_point() +
  geom_smooth(method = "lm") +
  xlab("Judgment Error") +
  ylab("EMD")

g3 = ggplot(df_to_show, mapping = aes(rt_err, look_dist)) +
  geom_point() +
  geom_smooth(method = "lm") + 
  xlab("RT Error") +
  ylab("EMD")

g = ggarrange(g1, g2, g3, ncol= 3)

`geom_smooth()` using formula 'y ~ x'
`geom_smooth()` using formula 'y ~ x'
`geom_smooth()` using formula 'y ~ x'

df_to_show = df_bandit_performance %>% 
  group_by(thresholds, 
           tradeoffs) %>% 
  summarise(ave_score = mean(combined_ranks))

ggplot(df_to_show, mapping = aes(x = thresholds,
                                 y = tradeoffs,
                                 fill = ave_score)) +
  geom_tile() +
  scale_y_continuous(trans="log10")

df_to_show = df_bandit_performance %>% 
  group_by(thresholds,
           bws) %>% 
  summarise(ave_score = mean(combined_ranks))

ggplot(df_to_show, mapping = aes(x = thresholds,
                                 y = bws,
                                 fill = ave_score)) +
  geom_tile()

df_to_show = df_bandit_performance %>% 
  group_by(thresholds,
           sample_weights) %>% 
  summarise(ave_score = mean(combined_ranks))

ggplot(df_to_show, mapping = aes(x = thresholds, 
                                 y = sample_weights, 
                                 fill = ave_score)) +
  geom_tile()

df_to_show = df_bandit_performance %>% 
  group_by(tradeoffs, bws) %>% 
  summarise(ave_score = mean(combined_ranks))

ggplot(df_to_show, mapping = aes(x = tradeoffs,
                                 y = bws,
                                 fill = ave_score)) +
  geom_tile() +
  scale_x_continuous(trans="log10")

df_to_show = df_bandit_performance %>% 
  group_by(tradeoffs, sample_weights) %>% 
  summarise(ave_score = mean(combined_ranks))

ggplot(df_to_show, mapping = aes(x = tradeoffs,
                                 y = sample_weights,
                                 fill = ave_score)) +
  geom_tile() +
  scale_x_continuous(trans = "log10")

df_to_show = df_bandit_performance %>% 
  group_by(bws, sample_weights) %>% 
  summarise(ave_score = mean(combined_ranks))

ggplot(df_to_show, mapping = aes(x = bws,
                                 y = sample_weights,
                                 fill = ave_score)) +
  geom_tile()

df_to_show = df_bandit_performance %>% 
  pivot_longer(c(thresholds, tradeoffs, bws, sample_weights),
               names_to = "param",
               values_to = "param_val") %>% 
  group_by(param, param_val) %>% 
  summarise(ave_score = mean(combined_ranks))# %>%
  # filter(!((param == "thresholds") & (param_val == 1.1)))

ggplot(df_to_show, mapping = aes(x = param_val, y = ave_score)) +
  geom_point() +
  geom_line() +
  facet_wrap(~param, scales = "free") +
  theme(axis.text.x = element_text(size = 10))

Sensitivity Analysis

df_bandit_performance %>% 
  filter(thresholds == 1.10)

# A tibble: 280 × 12
   thresholds tradeoffs   bws sample_weights judge_err rt_err     X look_dist
        <dbl>     <dbl> <dbl>          <dbl>     <dbl>  <dbl> <int>     <dbl>
 1        1.1     0.001    30           1000      43.6   17.6     0  0.000402
 2        1.1     0.001    40           1000      43.6   17.6     0  0.000402
 3        1.1     0.001    30           1050      43.6   17.6     0  0.000402
 4        1.1     0.001    40           1050      43.6   17.6     0  0.000402
 5        1.1     0.001    30           1100      43.6   17.6     0  0.000402
 6        1.1     0.001    40           1100      43.6   17.6     0  0.000402
 7        1.1     0.001    40            500      43.6   17.6     0  0.000402
 8        1.1     0.001    40            550      43.6   17.6     0  0.000402
 9        1.1     0.001    40            600      43.6   17.6     0  0.000402
10        1.1     0.001    30            650      43.6   17.6     0  0.000402
# … with 270 more rows, and 4 more variables: rank_judge <dbl>, rank_rt <dbl>,
#   rank_look_dist <dbl>, combined_ranks <dbl>

df_bandit_performance %>% 
  filter(thresholds != 1.10) %>% 
  group_by(thresholds) %>% 
  summarise(ave_judge_err = mean(judge_err),
            ave_rt_err = mean(rt_err),
            ave_look_dist = mean(look_dist)) %>% 
  ungroup() %>% 
  summarise(range_judge = range(ave_judge_err),
            range_rt = range(ave_rt_err),
            range_look_dist = range(ave_look_dist))

# max(temp) - min(temp)

# A tibble: 2 × 3
  range_judge range_rt range_look_dist
        <dbl>    <dbl>           <dbl>
1        21.1     9.72        0.000291
2        33.7    12.4         0.000320

df_bandit_performance %>% 
  filter(thresholds != 1.10) %>% 
  group_by(tradeoffs) %>% 
  summarise(ave_judge_err = mean(judge_err),
            ave_rt_err = mean(rt_err),
            ave_look_dist = mean(look_dist)) %>% 
  ungroup() %>% 
  summarise(range_judge = range(ave_judge_err),
            range_rt = range(ave_rt_err),
            range_look_dist = range(ave_look_dist))

# A tibble: 2 × 3
  range_judge range_rt range_look_dist
        <dbl>    <dbl>           <dbl>
1        23.2     10.3        0.000284
2        24.6     11.3        0.000323

df_bandit_performance %>% 
  filter(thresholds != 1.10) %>% 
  group_by(bws) %>% 
  summarise(ave_judge_err = mean(judge_err),
            ave_rt_err = mean(rt_err),
            ave_look_dist = mean(look_dist)) %>% 
  ungroup() %>% 
  summarise(range_judge = range(ave_judge_err),
            range_rt = range(ave_rt_err),
            range_look_dist = range(ave_look_dist))

# A tibble: 2 × 3
  range_judge range_rt range_look_dist
        <dbl>    <dbl>           <dbl>
1        22.2     10.2        0.000292
2        26.5     12.0        0.000328

df_bandit_performance %>% 
  filter(thresholds != 1.10) %>% 
  group_by(sample_weights) %>% 
  summarise(ave_judge_err = mean(judge_err),
            ave_rt_err = mean(rt_err),
            ave_look_dist = mean(look_dist)) %>% 
  ungroup() %>% 
  summarise(range_judge = range(ave_judge_err),
            range_rt = range(ave_rt_err),
            range_look_dist = range(ave_look_dist))

# A tibble: 2 × 3
  range_judge range_rt range_look_dist
        <dbl>    <dbl>           <dbl>
1        22.5     10.7        0.000304
2        24.3     10.8        0.000309

Fixed Sample

Load Heatmap performance

path = "../python/model/model_performance/grid_regression_error/"
reg_files = list.files(path)

df_reg_perf = tibble()

for (file in reg_files) {
  if (grepl("fixed_sample", file)) {
    temp = read.csv(paste(path, file, sep = ""))
    df_reg_perf = rbind(df_reg_perf, temp)
  }
}

df_reg_perf = df_reg_perf %>% 
  select(-X) %>% 
  rename(look_dist = sq_err)

num_samples_vec = c()
fixed_bws = c()
fixed_judge_err = c()
fixed_rt_err = c()
fixed_look_dist = c()

path = "../python/model/model_performance/grid_judgment_rt/"
list_filenames = list.files(path)


for (file in list_filenames) {
  file_list = str_split(file, "_")[[1]]
  if ("fixed" %in% file_list) {
    
    num_samples = as.numeric(file_list[match("samples", file_list) + 1])
    fixed_bw = as.numeric(file_list[match("bw", file_list) + 1])
    
    df_model_raw = read.csv(paste(path, file, sep = ""))
    
    judge_err = compare_fixed_judgment(df_data_mean_judge_train, df_model_raw)
    rt_err = compare_fixed_rt(df_data_mean_rt_train, df_model_raw)
    
    num_samples_vec = append(num_samples_vec, num_samples)
    fixed_bws = append(fixed_bws, fixed_bw)
    fixed_judge_err = append(fixed_judge_err, judge_err)
    fixed_rt_err = append(fixed_rt_err, rt_err)
    # fixed_look_dist = append(fixed_look_dist, dist)
  }
}

df_fixed_performance = tibble(num_samples = num_samples_vec,
                              bws = fixed_bws,
                              judge_err = fixed_judge_err,
                              rt_err = fixed_rt_err)

df_fixed_performance = df_fixed_performance %>% 
  left_join(df_reg_perf,
            by = c("num_samples",
                   "bws"))

df_fixed_performance = df_fixed_performance %>% 
  mutate(rank_judge = rank(judge_err),
         rank_rt = rank(rt_err),
         rank_look_dist = rank(look_dist),
         combined_rank = (rank_judge + rank_rt + rank_look_dist)/3) %>% 
  arrange(combined_rank)

df_to_show = df_fixed_performance

ggplot(data = df_to_show, mapping = aes(x = num_samples,
                                                  y = bws,
                                                  fill = combined_rank)) +
  geom_tile()

df_to_show = df_fixed_performance %>% 
  # filter(num_samples %in% seq(10,150,10),
  #        bws %in% seq(2,20,2)) %>%
  pivot_longer(c(num_samples, bws),
               names_to = "param",
               values_to = "param_val") %>% 
  group_by(param, param_val) %>% 
  summarise(ave_score = mean(combined_rank))


ggplot(df_to_show, mapping = aes(x = param_val,
                                 y = ave_score)) +
  geom_line() +
  geom_point() +
  facet_wrap(~param, scales = "free")

df_fixed_performance = df_fixed_performance %>%
  arrange(combined_rank)

Results Visualization

Bandit Model

df_model_judge_rt = read.csv("../python/model/model_performance/grid_judgment_rt/bandit_runs_30_threshold_0.95_tradeoff_0.003_sample_weight_950_bw_30.0_noise_params_0.2_0.8_0.2_trial_0_150.csv") %>% select(-X)

Judgments

df_model_mean_judge = df_model_judge_rt %>% 
  mutate(judgment = judgment + 1,
         judgment=factor(judgment)) %>%
  group_by(trial, judgment) %>% 
  summarise(model_mean = n()/(max(run)+1)) %>% 
  ungroup() %>% 
  complete(trial, judgment,
           fill = list(model_mean=0))

df_data_mean_judge_full = df_data_judge %>% 
  mutate(hole1 = as.numeric(response == 1),
         hole2 = as.numeric(response == 2),
         hole3 = as.numeric(response == 3)) %>% 
  select(-response) %>% 
  pivot_longer(c(hole1, hole2, hole3),
               names_to = "hole",
               values_to = "response") %>% 
  mutate(response = response) %>% 
  group_by(trial, hole) %>% 
  do(data.frame(rbind(smean.cl.boot(.$response)))) %>% 
  rename(human_mean = Mean,
         lower = Lower,
         upper = Upper)

df_human_mean_judge = df_data_mean_judge_full %>%
  mutate(hole = as.factor(str_sub(hole, -1, -1))) %>%
  rename(judgment = hole)

df_to_show = left_join(df_model_mean_judge,
                       df_human_mean_judge, 
                       by=c("trial", "judgment")) %>% 
  mutate(model = "Bandit") 

scaled_model = lm(human_mean ~ model_mean,
                  data = df_to_show)

scaled_model_predictions = predict(scaled_model)

model_cor = round(cor(scaled_model_predictions, df_to_show$human_mean), digits=2)
model_rmse = round(RMSE(scaled_model_predictions, df_to_show$human_mean), digits=2)

ggplot(data = df_to_show, mapping = aes(x = model_mean,
                                        y=human_mean)) +  
  geom_abline(slope = 1,
              intercept = 0,
              linetype="dotted") +
  geom_linerange(mapping = aes(ymin = lower,
                               ymax = upper),
                 alpha = 0.2) +
  geom_point(alpha=0.3) +
  geom_smooth(method = "lm",
              formula = y ~ x) +
  facet_grid(~ model) +
  xlab("Model Prediction") +
  ylab("Participant Response Proportion") +
  annotate("text",
           label = paste("r: ", model_cor),
           x=0.0,
           y=1,
           hjust=0) +
  annotate("text",
           label = paste("rmse: ", model_rmse),
           x=0.0,
           y=0.95,
           hjust = 0) +
  theme(plot.title = element_text(size=20, hjust=0.5),
        axis.title = element_text(size=16),
        axis.text = element_text(size=10))



# ggsave("figures/bandit_judgment.jpg", height = 4, width = 5)

temp = df_to_show$model_mean
sum((temp == 0) | (temp == 1))

[1] 358

# Number of responses on the end
temp = df_model_mean_judge$model_mean 
((temp == 1) ) %>% sum()

[1] 107

Response Time

df_model_mean_rt = df_model_judge_rt %>% 
  mutate(log_cols = ifelse(num_cols != 0, log(num_cols), num_cols)) %>% 
  group_by(trial) %>% 
  summarise(mean_time = mean(num_cols),
            mean_log_time = mean(log_cols))

df_data_mean_rt = df_data_rt %>%
  group_by(trial) %>% 
  summarise(mean_rt = mean(rt),
            mean_log_rt = mean(log(rt)))

df_to_show = left_join(df_model_mean_rt,
                       df_data_mean_rt,
                       by = c("trial"))

model_cor = round(cor(df_to_show$mean_time, df_to_show$mean_rt), digits=2)
model_rmse = round(RMSE(df_to_show$mean_time, df_to_show$mean_rt), digits=2)

ggplot(data = df_to_show, mapping = aes(x = mean_time, y = mean_rt)) +
  geom_point(alpha = 0.7,
             shape=16) +
  geom_smooth(method = "lm",
              formula = y ~ x) +
  # geom_label(mapping = aes(label = trial)) +
  ggtitle("Bandit Response Time") +
  xlab("Model Mean Collisions Across Runs") +
  ylab("Participant Mean Log Response Time") +
  annotate("text",
           label = paste("r =", model_cor),
           x=12,
           y=2500) +
  annotate("text",
           label = paste("rmse =", model_rmse),
           x=12,
           y=2200) +
  theme(plot.title = element_text(size=20, hjust=0.5),
        axis.title = element_text(size=14),
        axis.text = element_text(size=12))

# ggsave("figures/bandit_rt.png", height = 4, width = 5)

df_data_mean_rt = df_data_rt %>% 
  group_by(trial) %>% 
  do(data.frame(rbind(smean.cl.boot(.$log_rt)))) %>% 
  rename(mean_log_rt = Mean,
         upper = Upper,
         lower = Lower)

df_to_show = left_join(df_model_mean_rt,
                       df_data_mean_rt,
                       by = c("trial")) %>% 
  mutate(model = "Sequential Sampler")


# df_to_show = df_to_show %>% 
#   filter(mean_log_time > 0.5)


scaled_model = lm(mean_log_rt ~ 1 + mean_log_time,
                  data = df_to_show)

scaled_model_predictions = predict(scaled_model)

model_cor = round(cor(scaled_model_predictions, df_to_show$mean_log_rt), digits=2)
model_rmse = round(RMSE(scaled_model_predictions, df_to_show$mean_log_rt), digits=2)

ggplot(data = df_to_show, mapping = aes(x = mean_log_time, y = mean_log_rt)) +
  geom_linerange(mapping = aes(ymin = lower,
                               ymax = upper),
                 alpha = 0.15) +
  geom_point(alpha = 0.7,
             shape=16) +
  geom_smooth(method = "lm",
              formula = y ~ x) +
  facet_grid(~ model) +
  xlab("Model Mean Log Collisions") +
  ylab("Mean Log Response Time") +
  annotate("text",
           label = paste("r: ", model_cor),
           size = 6,
           x=0,
           y=8.8,
           hjust=0) +
  annotate("text",
           label = paste("rmse: ", model_rmse),
           size =6,
           x=0,
           y=8.68,
           hjust=0) +
  theme(plot.title = element_text(size=20, hjust=0.5),
        axis.title = element_text(size=20),
        axis.text = element_text(size=16),
        plot.margin = margin(10, 1, 1, 10))

ggsave("figures/bandit_log_rt.pdf", height = 4, width = 5)

df_bandit_rt = df_to_show %>% 
  rename(time_measure = mean_time,
         log_time = mean_log_time)

Fixed Sample Model

df_fixed_sample_judge_rt = read.csv("../python/model/model_performance/grid_judgment_rt/fixed_sample_num_samples_40_bw_50.0_noise_params_0.2_0.8_0.2_trial_0_150.csv") %>%  select(-X)

Judgments

df_fixed_sample_long = df_fixed_sample_judge_rt %>% 
  select(trial, hole1, hole2, hole3) %>% 
  pivot_longer(c(hole1, hole2, hole3),
               names_to = "hole",
               values_to = "prediction") 
  # mutate(judgment = factor(judgment))

df_to_show = df_fixed_sample_long %>% 
  left_join(df_data_mean_judge_full, by = c("trial", "hole")) %>% 
  mutate(model = "Fixed Sample")

fixed_sample_cor = round(cor(df_to_show$prediction, df_to_show$human_mean), digits = 2)
fixed_sample_rmse = round(RMSE(df_to_show$prediction, df_to_show$human_mean), digits = 2)

ggplot(df_to_show, mapping = aes(x = prediction, y = human_mean)) +
  geom_abline(slope = 1,
              intercept = 0,
              linetype = "dotted") +
  geom_linerange(mapping = aes(ymin = lower,
                               ymax = upper),
                 alpha=0.2) +
  geom_point(alpha=0.5,
             shape=16) +
  geom_smooth(method = "lm",
              formula = y ~ x) +
  annotate("text",
           label = paste("r:", fixed_sample_cor),
           x = 0.0,
           y = 1,
           hjust = 0) +
  annotate("text",
           label = paste("rmse:", fixed_sample_rmse),
           x = 0.0,
           y = 0.95,
           hjust = 0) +
  facet_grid(~ model) + 
  xlab("Model Prediction") +
  ylab("Participant Mean Judgment") +
  theme(plot.title = element_text(size=20, 
                                  hjust=0.5),
        axis.title = element_text(size=16),
        axis.text = element_text(size=10))

# ggsave("figures/fixed_sample_judgments.jpg", height=4, width=5)

temp = df_to_show$prediction %>% round(digits = 6)
sum((temp == 0) | (temp == 1))

[1] 187

df_fixed_judge = df_to_show %>% 
  mutate(judgment = as.factor(str_sub(hole, -1, -1)),
         model = "Uniform Sampler") %>% 
  select(-hole)

Response Time

df_to_show = df_fixed_sample_judge_rt %>% 
  select(trial, num_cols) %>% 
  mutate(time_measure = num_cols,
         log_time = log(num_cols)) %>% 
  left_join(df_data_mean_rt, by = "trial") %>% 
  mutate(model = "Uniform Sampler")


scaled_model = lm(mean_log_rt ~ 1 + log_time,
                  data = df_to_show)

scaled_model_predictions = predict(scaled_model)

fixed_sample_rt_cor = round(cor(scaled_model_predictions, df_to_show$mean_log_rt), digits = 2)
fixed_sample_rt_rmse = round(RMSE(scaled_model_predictions, df_to_show$mean_log_rt), digits = 2)

ggplot(data = df_to_show, mapping = aes(x = log_time, y = mean_log_rt)) +
  geom_linerange(mapping = aes(ymin = lower,
                               ymax = upper),
                 alpha = 0.15) +
  geom_point(alpha=0.7,
             shape=16) +
  geom_smooth(method = "lm",
              formula = y ~ x) +
  facet_grid(~ model) +
  xlab("Model Log Collisions") +
  ylab("Mean log Response Time") +
  annotate("text",
           label = paste("r:", fixed_sample_rt_cor),
           x=5.1,
           y=8.8,
           size=6,
           hjust=0) +
  annotate("text",
           label = paste("rmse:", fixed_sample_rt_rmse),
           x=5.1,
           y=8.68,
           size=6,
           hjust=0) +
  theme(plot.title = element_text(size=20, hjust=0.5),
        axis.title = element_text(size=20),
        axis.text = element_text(size=16),
        plot.margin = margin(10,1,1,10))

ggsave("figures/fixed_sample_rt.pdf", height=4, width = 5)

Cogsci Figures

Judgments

df_bandit_judge = df_model_mean_judge %>% 
  mutate(judgment = paste("hole", judgment, sep = ""),
         model = "Sequential Sampler") %>% 
  rename(prediction = model_mean,
         hole = judgment)

df_fixed_judge = df_fixed_sample_long %>% 
  mutate(model = "Uniform Sampler")

df_to_show = rbind(df_bandit_judge,
                   df_fixed_judge) %>% 
  left_join(df_data_mean_judge_full)

Joining, by = c("trial", "hole")

df_sum_stat = df_to_show %>% 
  group_by(model) %>% 
  summarise(r = round(cor(prediction, human_mean), digits = 2),
            rmse = round(RMSE(prediction, human_mean), digits = 2))

ggplot(df_to_show, mapping = aes(x = prediction, y = human_mean)) +
  geom_abline(slope = 1,
              intercept = 0,
              linetype = "dotted") +
  geom_linerange(mapping = aes(ymin = lower,
                               ymax = upper),
                 alpha = 0.2) +
  geom_point(alpha = 0.2) +
  geom_smooth(method = "lm") + 
  geom_text(data = df_sum_stat,
            x = 0.0,
            y = 1,
            size = 6,
            hjust = 0,
            mapping = aes(label = paste("r: ", r, sep = ""))) +
  geom_text(data = df_sum_stat,
            x = 0.0,
            y = 0.93,
            size = 6,
            hjust = 0,
            mapping = aes(label = paste("rmse: ", rmse, sep = ""))) +
  facet_wrap(~ model) +
  scale_x_continuous(breaks = c(0.0, 0.25, 0.50, 0.75, 1.00),
                     labels = c("0%", "25%", "50%", "75%", "100%")) +
  scale_y_continuous(breaks = c(0.0, 0.25, 0.50, 0.75, 1.00),
                     labels = c("0%", "25%", "50%", "75%", "100%")) +
  xlab("Model Prediction") + 
  ylab("Participant Selection") +
  theme(plot.title = element_text(size=20, hjust=0.5),
        axis.title = element_text(size=24),
        axis.text = element_text(size=16),
        panel.spacing = unit(2, "lines"))

`geom_smooth()` using formula 'y ~ x'

ggsave("figures/model_judgment.pdf",
       width = 10,
       height = 4)

`geom_smooth()` using formula 'y ~ x'

cor(df_bandit_judge$prediction, df_fixed_judge$prediction)

[1] 0.899042

EMD

df_emd_bandit = read.csv("../python/model/model_performance/emd/bandit.csv") %>% 
  select(trial, distance) %>% 
  mutate(model = "Sequential Sampler")

df_emd_fixed_sample = read.csv("../python/model/model_performance/emd/fixed_sample.csv") %>%
  select(trial, distance) %>%
  mutate(trial = factor(trial),
         model = "Uniform Sampler")

df_emd_baseline = read.csv("../python/model/model_performance/emd/visual_features.csv") %>% 
  select(-X) %>% 
  mutate(trial = factor(trial),
         model = "Visual Features")

to_highlight = c()

set.seed(1)

df_to_show = rbind(df_emd_bandit, df_emd_fixed_sample, df_emd_baseline) %>% 
  mutate(model = factor(model,
                        levels = c("Sequential Sampler", "Uniform Sampler", "Visual Features"),
                        labels = c(1,2,3)),
         model = as.numeric(as.character(model)),
         highlight = trial %in% to_highlight,
         model_jitter = model + runif(n = n(),
                                      min = -0.15,
                                      max = 0.15)) 

# ggplot(df_to_show, mapping = aes(x = model, y = distance)) +
ggplot(df_to_show, mapping = aes(x = model, 
                                 y = distance, 
                                 color = highlight)) +
  geom_line(mapping = aes(x = model_jitter, group = trial), 
            alpha = 0.05) +
  geom_point(mapping = aes(x = model_jitter),
             alpha = 0.5,
             shape=16,
             size=3) +
  stat_summary(fun.data = "mean_cl_boot", color = "red", size=0.8) +
  scale_x_continuous(breaks = c(1,2,3), labels = c("Sequential Sampler", "Uniform Sampler", "Visual Features")) + 
  scale_color_manual(values = c("black", "magenta3")) +
  ylab("Earth Mover's Distance") +
  theme(legend.title = element_blank(),
        legend.position = "none",
        axis.title.y = element_text(size=24),
        axis.title.x = element_blank(),
        axis.text = element_text(size=16))

ggsave("figures/emd_comparison.pdf",
       height = 5,
       width = 8)

df_emd = rbind(df_emd_bandit,
               df_emd_fixed_sample,
               df_emd_baseline)

df_emd %>% 
  group_by(model) %>% 
  do(data.frame(rbind(round(smean.cl.boot(.$distance), 2))))

# A tibble: 3 × 4
# Groups:   model [3]
  model               Mean Lower Upper
  <chr>              <dbl> <dbl> <dbl>
1 Sequential Sampler  49.2  46.9  51.5
2 Uniform Sampler     65.0  62.6  67.6
3 Visual Features     71.5  69.0  73.9